assumptions
Humboldt-Universität zu Berlin
2023-04-13
Today we will learn…
Sections 6.3, 6.4 and 7.9 in Winter (2019)
Start with a clean R Environment (Session > Restart R).
Suppress scientific notation (make very small numbers easier to read):
# load in dataset
df_crit_verb <-
readr::read_csv(
here::here("data/tidy_data_lifetime_pilot.csv"),
# for special characters
locale = readr::locale(encoding = "latin1")
) |>
mutate_if(is.character, as.factor) |> # all character variables as factor
mutate(lifetime = fct_relevel(lifetime, "living", "dead"),
tense = fct_relevel(tense, "PP", "SF")) |>
filter(type == "critical", # only critical trials
px != "px3", # px3 had a lot of missing values
fp > 0, # only values of fp above 0
region == "verb") %>% # critical region only
droplevels() # remove any factor levels with no observations\[ \begin{align} observed\;values &= fitted\;values + residuals\\ residuals &= observed\;values - fitted\;values \end{align} \]
Set contrasts
Check contrasts
Call:
lm(formula = fp ~ lifetime * tense, data = .)
Residuals:
Min 1Q Median 3Q Max
-240.62 -108.69 -27.62 56.20 778.65
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 309.06 6.26 49.367 <0.0000000000000002 ***
lifetime1 31.52 12.52 2.517 0.0121 *
tense1 -12.75 12.52 -1.018 0.3090
lifetime1:tense1 -21.69 25.04 -0.866 0.3868
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 145.9 on 539 degrees of freedom
Multiple R-squared: 0.01499, Adjusted R-squared: 0.009506
F-statistic: 2.734 on 3 and 539 DF, p-value: 0.04306
Image source: Winter (2019) (all rights reserved)
Figure 1: Visualising normality of residuals
for more see Section 5.4 in Winter (2019)
the R funtion log() computes the ‘natural logarithm’ (and is the inverse of the exponential exp())
log() makes large numbers smaller
exp() makes small numbers larger
[1] 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379 1.7917595 1.9459101
[8] 2.0794415 2.1972246 2.3025851
[1] 2.718282 7.389056 20.085537 54.598150 148.413159
[6] 403.428793 1096.633158 2980.957987 8103.083928 22026.465795
tt < 500ms), with some larger values (> tt 1000)
Call:
lm(formula = log(fp) ~ lifetime * tense, data = .)
Residuals:
Min 1Q Median 3Q Max
-1.18141 -0.34282 0.00058 0.26923 1.38822
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.63477 0.01877 300.193 < 0.0000000000000002 ***
lifetime1 0.10130 0.03754 2.698 0.00718 **
tense1 -0.03357 0.03754 -0.894 0.37158
lifetime1:tense1 -0.08320 0.07508 -1.108 0.26828
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4374 on 539 degrees of freedom
Multiple R-squared: 0.01712, Adjusted R-squared: 0.01165
F-statistic: 3.129 on 3 and 539 DF, p-value: 0.02538
Raw first-pass reading times
augment(fit_fp_log, data = df_crit_verb[df_crit_verb$fp > 0,]) %>%
distinct(lifetime,tense,.fitted) %>%
arrange(lifetime) # A tibble: 4 × 3
lifetime tense .fitted
<fct> <fct> <dbl>
1 living PP 5.58
2 living SF 5.59
3 dead PP 5.72
4 dead SF 5.65
df_crit_verb %>%
mutate(predicted_raw = predict(fit_fp),
predicted_log = predict(fit_fp_log)) %>%
distinct(lifetime,tense,predicted_raw, predicted_log) %>%
arrange(lifetime) # A tibble: 4 × 4
lifetime tense predicted_raw predicted_log
<fct> <fct> <dbl> <dbl>
1 living PP 294. 5.58
2 living SF 292. 5.59
3 dead PP 337. 5.72
4 dead SF 313. 5.65
(Intercept)
279.9937
(Intercept)
309.0606
(Intercept)
293.302
(Intercept)
266.1646
(Intercept)
294.5413
# living-SF
exp(coef(fit_fp_log)['(Intercept)']) * exp(coef(fit_fp_log)['lifetime1'] * -0.5 + coef(fit_fp_log)['tense1'] * -0.5 )(Intercept)
270.6701
lifetime or tense, and would expect them to replicate if we were to re-run the studylifetime, tense), and random effects (in eye-tracking experiments, typically participant, item)
performance packageperformance package has a nice function check_model() that produces 6 plots to check model fitFigure 2: Performance of model with log reading times
Figure 3: Performance of model with log reading times
tidy() function from broom is handy in combination with knitr::kable(digits = x) and kableExtra::kable_styling()
tidy(fit_fp_log) %>%
kable(digits = 10,
col.names = c("Coefficient",
"Estimate (log)",
"SE",
"t-value",
"p-value")) %>%
kable_styling()| Coefficient | Estimate (log) | SE | t-value | p-value |
|---|---|---|---|---|
| (Intercept) | 5.63476711 | 0.01877047 | 300.1931843 | 0.000000000 |
| lifetime1 | 0.10130415 | 0.03754094 | 2.6984981 | 0.007183743 |
| tense1 | -0.03357132 | 0.03754094 | -0.8942589 | 0.371582493 |
| lifetime1:tense1 | -0.08320391 | 0.07508188 | -1.1081757 | 0.268280198 |
# source: https://stackoverflow.com/questions/37470202/in-line-code-for-reporting-p-values-001-in-r-markdown
# OR USE
pacman::p_load(broman)
format_pval <- function(x){
if (x < .001) return(paste('<', '.001'))
if (x < .01) return(paste('<', '.01'))
if (x < .05) return(paste('<', '.05'))
paste('=', myround(x, 3)) # if above .05, print p-value to 3 decimalp points
}tidy(fit_fp_log) %>%
mutate(format = sapply(p.value, function(x) format_pval(x))) %>%
mutate(signif = sapply(p.value, function(x) make_stars(x))) %>%
select(-p.value) %>%
kable(digits = 10,
col.names = c("Coefficient",
"Estimate (log)",
"SE",
"t-value",
"p-value",
"sign")) %>%
kable_styling()| Coefficient | Estimate (log) | SE | t-value | p-value | sign |
|---|---|---|---|---|---|
| (Intercept) | 5.63476711 | 0.01877047 | 300.1931843 | < .001 | *** |
| lifetime1 | 0.10130415 | 0.03754094 | 2.6984981 | < .01 | ** |
| tense1 | -0.03357132 | 0.03754094 | -0.8942589 | = 0.372 | |
| lifetime1:tense1 | -0.08320391 | 0.07508188 | -1.1081757 | = 0.268 |
Figure 4: ?(caption)
Figure 5: ?(caption)
fig_error <-
df_crit_verb %>%
filter(fp > 0) %>%
summarySEwithin(measurevar="fp", withinvars=c("lifetime", "tense"), idvar="px") %>%
mutate(upper = fp+ci,
lower = fp-ci) %>%
ggplot(aes(x = lifetime, y = fp, colour = tense, shape = tense)) +
labs(title = "Mean first-pass reading times (with 95% CIs)") +
geom_point(position = position_dodge(0.2), size = 2) +
geom_line(position = position_dodge(0.2), aes(group=tense)) +
geom_errorbar(aes(ymin=lower,ymax=upper), position = position_dodge(0.2), width = .2) +
theme_bw() +
theme(text = element_text(size=8))
fig_error_log <-
df_crit_verb %>%
filter(fp > 0) %>%
mutate(fp_log = log(fp)) %>%
summarySEwithin(measurevar="fp_log", withinvars=c("lifetime", "tense"), idvar="px") %>%
mutate(upper = fp_log+ci,
lower = fp_log-ci) %>%
ggplot(aes(x = lifetime, y = fp_log, colour = tense, shape = tense)) +
labs(title = "Mean first-pass reading times (with 95% CIs)") +
geom_point(position = position_dodge(0.2), size = 2) +
geom_line(position = position_dodge(0.2), aes(group=tense)) +
geom_errorbar(aes(ymin=lower,ymax=upper), position = position_dodge(0.2), width = .2) +
theme_bw() +
theme(text = element_text(size=8))Figure 6: ?(caption)
we looked at some assumptions of linear models and how to check them
we saw how to log-transform our data so that our residuals are normally distributed
we learned the importance of the independence assumption, which is violated by repeated-measures designs
next, we’ll learn how to account for lack of independence between the data with mixed models
Model assumptions